Code
sshfs narecco@ant-login.linux.crg.es:/users/mirimia ~/mntLast code execution: 2024 February 07, Wednesday @ 19:02:17.
We performed histone PTMs DIA mass spectrometry and I analysed the raw files with EpiProfile v2.1 (Yuan et al. 2018). Here I parse the EpiProfile output and analyse the PTMs values with DEP (Zhang et al. 2018).
Acknowledgements: I’d like to thank Zuo-Fei Yuan and Simone Sidoli from the Ben Garcia lab for their help with data analysis and sample preparation respectively. Also a big thank you to Cristina Chiva from the CRG/UPF proteomics facility.
To fetch files on the CRG cluster I mount the server on my local computer with the following command in the terminal.
sshfs narecco@ant-login.linux.crg.es:/users/mirimia ~/mntLoad packages required for the analysis and suppress any message. Check the Section 6 section at the end for package versions.
library(dplyr, warn.conflicts = F, quietly = T)
library(readr)
library(readxl)
library(pheatmap)
library(DEP)
library(Cairo)
library(DT)Load my R package niar to use my custom made function to parse EpiProfile output.
if ( ('niar' %in% .packages(all.available = TRUE)) == TRUE ) {
library('niar')
} else if ( ('niar' %in% .packages(all.available = TRUE)) == FALSE ){
message("The package niar is not installed so I'm trying to load it from ",
"the local repository of the package")
if ( dir.exists('~/mnt/narecco/software/R/niar' ) ) {
devtools::load_all(path = '~/mnt/narecco/software/R/niar')
} else {
stop("Can't find the local repo of the niar package! ",
"You must install it with:\n",
"devtools::install_github('Ni-Ar/niar') ")
}
} else{
stop("Can't understand if 'niar' package was installed beforehand")
}The package niar is not installed so I'm trying to load it from the local repository of the package
ℹ Loading niar
Define a plot style.
theme_classic(base_size = 6.3, base_family = "Arial") +
theme(plot.title = element_text(size = 8, vjust = -1, hjust = 0),
plot.background = element_blank(),
strip.background = element_blank(),
panel.background = element_blank(),
axis.line = element_line(linewidth = 0.1),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 5),
axis.text = element_text(colour = "black"),
axis.text.x = element_text(margin = margin(t = 0, unit = "mm")),
axis.ticks.x = element_blank(),
axis.ticks.y = element_line(linewidth = 0.1),
axis.ticks.length.y = unit(1, "mm"),
legend.margin = margin(l = 0, t = 0, unit = "mm"),
legend.box.margin = margin(l = 0, unit = "mm"),
legend.box.background = element_blank(),
legend.box.spacing = unit(0.5, "mm"),
# legend.title = element_blank(),
legend.text = element_text(margin = margin(l = -10, unit = "mm")),
legend.text.align = 0,
legend.key.height = unit(2, "mm"),
legend.key.width = unit(7.5, "mm"),
legend.key = element_blank(),
legend.background = element_blank(),
legend.spacing.x = unit(3, 'mm')
) -> alluvial_thm I modify the geom_stratum() function just to add a very small line width the the black borders of the columns. With linewidth = .05, I don’t have to change them manually later in Illustrator. This will hopefully fixed in future
Small function to parse the PTM abundance as 10 to the power scientific notion.
scientific_10 <- function(x) {
parse(text=gsub("e\\+", " %*%10^", scales::scientific_format()(x)))
}Here I organise all the file and directories paths I need to run the analysis and define where to save the processed tables and figures.
oneDrive_Dir <- file.path("~/OneDrive - CRG - Centre de Regulacio Genomica/Suz12_AS_project")
code_dir_fig <- file.path(oneDrive_Dir, "_Code/Fig4")
tbl_dir_fig4 <- file.path(code_dir_fig, "tables")
pdf_dir_fig4 <- file.path(code_dir_fig, "pdfs")
if (!dir.exists(pdf_dir_fig4)) { dir.create(pdf_dir_fig4, recursive = T) }
if (!dir.exists(tbl_dir_fig4)) { dir.create(tbl_dir_fig4, recursive = T) }Path to medata data and EpiProfile output files.
metadata_path <- file.path(tbl_dir_fig4, 'Histone_PTMs_metadata.xlsx')
WT_Dex4_histone_ratio_path <- list.files(tbl_dir_fig4,
pattern = "FullSet1_Run2_DIA_histone_ratios_Suz12dExon4.xls",
full.names = T)
Resce_KO_histone_ratio_path <- list.files(tbl_dir_fig4,
pattern = "FullSet2_Run1_DIA_histone_ratios_Suz12KO.xls",
full.names = T)
stopifnot(file.exists(metadata_path))
stopifnot(file.exists(WT_Dex4_histone_ratio_path))
stopifnot(file.exists(Resce_KO_histone_ratio_path))SUZ12 WT proteoforms quantification data.
protein_psi_path <- file.path(tbl_dir_fig4, 'Protein_PSI_3reps_LS.xlsx')Import metadata.
mtdt <- read_excel(metadata_path) |>
setNames(c("Sample", "Short_Name", "Clone_Name", "Condition", "Replicate"))
mtdt$Condition <- factor(mtdt$Condition,
levels = c("WT", "CSexon4", "Dexon4", "KO", "KO_L", "KO_S") )Import DIA data quantification using automatic EpiProfile v2.1 quantifications.
read_EpiProfile_histone_ratio(WT_Dex4_histone_ratio_path) |>
left_join(y = mtdt, by = "Sample") |>
mutate(Facility_ID = Sample) |>
mutate(Sample = Short_Name) |>
tidy_hPTMs(split_double_PTMs = T) -> dia_auto_set1
read_EpiProfile_histone_ratio(Resce_KO_histone_ratio_path) |>
# fix an extra _ in the file name
mutate(Sample = gsub(pattern = '2022MH003_NIAR_005_02__30pto_DIA',
replacement = '2022MH003_NIAR_005_02_30pto_DIA',
x = Sample) ) |>
left_join(y = mtdt, by = "Sample") |>
mutate(Facility_ID = Sample) |>
mutate(Sample = Short_Name) |>
tidy_hPTMs(split_double_PTMs = T) -> dia_auto_set2
data_dia_auto <- rbind(dia_auto_set1, dia_auto_set2)Tidy up the information for histone peptide H3 K27- K36. Here I split the H3 K27 and K36 signal coming form the same peptide. In this way I can pile up the all methyl-forms.
data_dia_auto |>
subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
mutate(K27 = case_when( is.na(Sub_PTM1) & grepl(pattern = "^K27", x = Modification) ~ Modification,
grepl(pattern = "unmod", x = Modification) ~ "unmod" ) ) |>
mutate(K36 = case_when( is.na(Sub_PTM2) & grepl(pattern = "^K36", x = Modification) ~ Modification,
grepl(pattern = "unmod", x = Modification) ~ "unmod" ) ) |>
relocate(K27, K36, .after = Sub_PTM2) |>
unite(col = "K27", c("Sub_PTM1", "K27"), sep = "", na.rm = T, remove = T ) |>
unite(col = "K36", c("Sub_PTM2", "K36"), sep = "", na.rm = T, remove = T ) |>
mutate(K27 = case_when(grepl(pattern = "^$", x = K27) ~ "unmod",
grepl(pattern = "[aA-zZ]", x = K27) ~ K27 )) |>
mutate(K36 = case_when(grepl(pattern = "^$", x = K36) ~ "unmod",
grepl(pattern = "[aA-zZ]", x = K36) ~ K36 )) |>
# When one mark is not found I label it unmodified ('unmod') but it could also
# be labelled as "only the other mark".
# mutate(K27 = case_when(grepl(pattern = "^$", x = K27) ~ "Only\nK36me",
# grepl(pattern = "[aA-zZ]", x = K27) ~ K27 )) |>
# mutate(K36 = case_when(grepl(pattern = "^$", x = K36) ~ "Only\nK27ac-me",
# grepl(pattern = "[aA-zZ]", x = K36) ~ K36 )) |>
relocate(K27, K36, .after = PTM) |>
mutate(Histone = factor(Histone, levels = c("H3", "H3.3") ) ) -> data_K27_36Overview of H3 and H3.3 PTMs in lysines 27 and 36
datatable(data_K27_36, rownames = FALSE, filter = 'top',
options = list(pageLength = 7, autoWidth = TRUE) )I can also deconvolve the K27 and K36 signals from the same peptide and pile-up each signal.
data_K27_36 |>
group_by(Condition, K27, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) )|>
select(Histone, Condition, Average_Ratio, K27, Short_Name, Ratio) |>
mutate(Ratio_Sum = sum(Average_Ratio) / 3) |>
mutate(Percent = round(Average_Ratio*100, 2) ) |>
select(Condition, Histone, K27, Average_Ratio, Ratio_Sum, Percent) |>
unique() |>
mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
arrange(Condition, Histone, K27) |> ungroup() -> data_K27Same as above but for K36.
data_K27_36 |>
group_by(Condition, K36, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
mutate(K36 = factor(K36, levels = c("K36me3", "K36me2", "K36me1", "unmod") ) ) |>
select(Histone, PTM, Condition, Average_Ratio, K36, Short_Name, Ratio) |>
mutate(Ratio_Sum = sum(Average_Ratio) / 3) |>
mutate(Percent = round(Average_Ratio*100, 2) ) |>
select(Condition, Histone, K36, Average_Ratio, Ratio_Sum, Percent) |>
unique() |>
mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
arrange(Condition, Histone, K36) |> ungroup() -> data_K36For specific histone PTMs one can explore the results in the following tables.
Overview of H3 and H3.3 K27 acetylation and methylation quantifications
datatable(data_K27, rownames = FALSE, filter = 'top',
options = list(pageLength = 5, autoWidth = TRUE) )Overview of H3 and H3.3 K36 methylation quantifications
datatable(data_K36, rownames = FALSE, filter = 'top',
options = list(pageLength = 5, autoWidth = TRUE) )Plot figure 4B panel with the pileup of figure H3K27 methylation and acetylation.
data_K27_36 |>
group_by(Condition, K27, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3) |> # 3 replicates
subset(Histone %in% "H3") |>
mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) ) |>
ggplot() +
aes(x = Condition, stratum = K27, alluvium = K27,
y = Ratio_Sum, fill = K27) +
geom_alluvium(alpha = 0.5, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, colour = "black") +
geom_fit_text(aes(label = K27 ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = c("K27ac" = "#feac81", "K27me1" ="#ced1af",
"K27me2" = "#748f46", "K27me3" = "#47632a",
"unmod" = "gray" ),
name = 'H3') +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm -> p_Alluvial_K27
p_Alluvial_K27 Save to pdf.
ggsave(filename = "Fig4B_H3_only_K27_Mrgd_Alluvial.pdf", path = pdf_dir_fig4,
plot = p_Alluvial_K27, device = cairo_pdf, units = "cm",
width = 6.5, height = 4.85)Use the area under the peak as counts to be analysed using DEP framework.
rbind( read_EpiProfile_histone_ratio(WT_Dex4_histone_ratio_path),
read_EpiProfile_histone_ratio(Resce_KO_histone_ratio_path) ) |>
# fix an extra _ in the file name like before
mutate(Sample = gsub(pattern = '2022MH003_NIAR_005_02__30pto_DIA',
replacement = '2022MH003_NIAR_005_02_30pto_DIA',
x = Sample) ) |>
left_join(y = mtdt, by = "Sample") |>
mutate(Facility_ID = Sample) |>
mutate(Sample = Short_Name) -> data_dia_auto_not_tidyOne small caveat in this analysis is that the areas under the peaks are very large numbers or at least larger than 2147483647 which is the highest number R can store as integer. To convert the areas to integers from numeric type one could use the package bit64 that stores integers in 64 bits instead of the default 32 bit. My solution was however simpler, I divided the all values by 1000 to scale down by 3 order of magnitude. I then convert the dataframe to matrix.
I tried to add a pseudocounts to the data but it actually doesn’t help with the analysis and is better to input the areas as they are. The function make_unique is the first step of the analysis of DEP and creates a unique data frame.
data_dia_auto_not_tidy |>
tidy_hPTMs() |>
select(First_Col, PTM, Sample, Area) |>
# reduce all areas by a 1000 to avoid having to coerc to integer big numbers
mutate(Area = Area / 1e4) |>
pivot_wider(id_cols = c(First_Col, PTM), names_from = Sample, values_from = Area) |>
# Coerce the areas to integers
mutate(across(.cols = ends_with(c("_1", "_2", "_3")), .fns = as.integer )) |>
make_unique(names = "First_Col", ids = "PTM") -> data_uniqueThe data frame data_unique contains 22 columns (of which 18 are samples) and 185 rows corresponding to the identified histone PTMs.
Create a metadata compatible for DEP set the KO as the reference sample as I’m interesting in the rescue efficiency relative to the KO.
# Set metadata colnames to be as `make_se` wants them to be
colnames(mtdt) <- c("Sample", "label", "Clone_Name", "condition", "replicate")
mutate(mtdt, condition = factor(condition, levels = c("WT", "CSexon4", "Dexon4",
"KO", "KO_L", "KO_S"))) |>
print(n = 3) -> mtdt_dep# A tibble: 18 × 5
Sample label Clone_Name condition replicate
<chr> <chr> <chr> <fct> <dbl>
1 2022MZ006_IVMO_001_02_25pto_DIA WT_1 ZWT WT 1
2 2022MZ006_IVMO_002_02_25pto_DIA WT_2 ZA1 WT 2
3 2022MZ006_IVMO_003_02_25pto_DIA WT_3 ZA2 WT 3
# ℹ 15 more rows
Find the columns that contain the samples areas in the data_unique dataframe and create a SummarizedExperiment object to be used in for the DEP analysis. The first step is to do a log2 transformation of the area counts.
Samples_area_columns <- grep(pattern = paste(mtdt_dep$label, collapse = "|"), colnames(data_unique))
data_se <- make_se(data_unique, columns = Samples_area_columns, expdesign = mtdt_dep)Now the data_se contains the modifications in the rows (assuming they are proteins) and in the column the area.
data_filt <- filter_missval(data_se, thr = 0)Check PTMs frequency across samples.
plot_frequency(data_filt)Overall most modifications are found in all samples.
data_norm <- normalize_vsn(data_filt)Verify the fit.
meanSdPlot(x = data_norm)Seems okay.
Check normalisation effect on data.
plot_normalization(data_se, data_filt, data_norm) +
scale_fill_manual(values = c("WT" = "royalblue3", "Dexon4" = "goldenrod",
"CSexon4" = "#B0461C", "KO" = "gray45",
"KO_L" = "mediumpurple3", "KO_S" = "darkorange1")) +
theme(axis.text = element_text(color = 'black'),
plot.background = element_blank(),
panel.background = element_blank())Check missing values.
plot_missval(data_norm)I don’t see obvious missing values trends, besides only that replicate #3 of CSexon has less PTMs detected. I don’t think it’s necessary to impute missing values.
Test every sample versus WT control
data_diff <- test_diff(data_norm, type = "control", control = "WT")Significance thresholds:
P-value <= 0.05
|Log2 Fold Change| >= 1.5
dep <- add_rejections(data_diff, alpha = 0.05, lfc = 1.5 )Check sample clustering a Gower’s distance matrix.s
plot_dist(dep, significant = TRUE, pal = "PuOr") This is basically showing the 2 batches of the 2 MS acquisition experiments.
Check sample clustering using the same distance.
plot_heatmap(dep, type = "centered", kmeans = F, clustering_distance = "gower",
col_limit = 4, show_row_names = TRUE,
indicate = c("condition", "replicate"))Get results to make a volcano plot.
data_results <- get_results(dep)
data_results |>
as_tibble() |>
select(-significant) |>
pivot_longer(cols = !c("name", "ID"),
names_to = c("contrast", "Feature"),
names_sep = c("_(vs_WT_|centered)"),
values_to = "Value"
) |>
mutate(Feature = ifelse(test = Feature == "", yes = "centered", no = Feature)) |>
pivot_wider(names_from = "Feature", values_from = "Value") |>
dplyr::select(-ID) -> res_all_tidy
colnames(res_all_tidy)[colnames(res_all_tidy) == "name"] <- 'First_Col'Extract WT sample area of each PTM, to be used for the volcano plot.
data_dia_auto |>
subset(Condition == "WT") |>
group_by(First_Col) |>
mutate(WT_Mean_Area = mean(Area) + 1 ) |>
relocate(WT_Mean_Area, .after = PTM) |>
select(c(First_Col, WT_Mean_Area)) |>
unique() -> WT_areaGet the histone PTMs info columns.
histone_ptms_anno <- unique(data_dia_auto[, 1:10])Add all the info the the tidy results dataframe.
res_all <- left_join(x = res_all_tidy, y = histone_ptms_anno, by = join_by(First_Col))
res_all <- left_join(res_all, WT_area, by = join_by(First_Col) )
res_all |>
mutate(Histone = factor(Histone, levels = sort(unique(res_all$Histone)) ) ) |>
arrange(desc(Histone), Peptide_Start, Peptide_End) |>
ungroup() -> res_allExport all tidy results for supplementary table 3 as tab delimited, or as csv with UTF-8 Byte order mark which indicates to Excel the csv is UTF-8 encoded.
write_delim(x = res_all, file = file.path(tbl_dir_fig4, 'DEP_analysis_histone_MS.tab'),
delim = '\t', col_names = T, append = F, quote = 'none', na = "NA",
progress = F, escape = 'none')
write_excel_csv(x = res_all,
file = file.path(tbl_dir_fig4, 'DEP_analysis_histone_MS.csv'),
delim = ',', col_names = T, append = F, quote = 'none',
na = "NA", progress = F, escape = 'none')Interactively explore the table
res_all |>
mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
datatable(rownames = FALSE, filter = 'top',
options = list(pageLength = 3, autoWidth = TRUE) )Plot centred abundance
res_all |>
subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
subset(!is.na(centered)) |>
ggplot(aes(x = centered, y = PTM, fill = contrast)) +
facet_grid(rows = vars(Modification), scales = 'free' ) +
geom_vline(xintercept = 0, linewidth = 0.4) +
geom_col(width = 0.75,
position = position_dodge(width = 0.75, preserve = 'single')) +
geom_point(aes(size = WT_Mean_Area),
shape = 21, position = position_dodge(width = 0.75),
stroke = 0.2) +
labs(x = "Centered PTM abundance") +
scale_x_continuous(limits = c(-5, 5), oob = scales::oob_squish ) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
scale_fill_manual(values = c("WT" = "royalblue3", "Dexon4" = "goldenrod",
"CSexon4" = "#B0461C", "KO" = "gray45",
"KO_L" = "mediumpurple3", "KO_S" = "darkorange1"),
name = "Sample") +
theme_classic() +
theme(strip.text = element_blank(),
strip.background = element_blank(),
axis.text = element_text(colour = "black"),
axis.title.y = element_blank(),
legend.key.size = unit(1, 'mm'),
panel.grid.major.y = element_line()) -> p_centered_K27_K36
p_centered_K27_K36Show differential histone PTMs as volcano plots.
Prepare data for plot.
res_Dex4 <- subset(res_all, contrast == "Dexon4" & !is.na(p.val))
res_Dex4 <- res_Dex4 |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) Plot figure 4C.
ggplot(res_Dex4 ) +
aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
geom_point(shape = 21, stroke = 0.2) +
annotate(geom = "label", x = 5.5, y = 0.10, label = "∆ex4",
colour = "black", fill = "goldenrod", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
annotate(geom = "label", x = -5.5, y = 0.10, label = "WT",
colour = "white", fill = "royalblue3", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
# points on the right side
geom_label_repel(data = subset(res_Dex4, p.adj < 0.09 & ratio > 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = -0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
# points on the left side
geom_label_repel(data = subset(res_Dex4, p.adj < 0.1 & ratio < 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = 0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
guides(size = guide_legend(override.aes = list(fill = "white"))) +
scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-6.4, 6.4), n.breaks = 7) +
scale_y_continuous(expand = expansion(add = c(0, 0.25), mult = 0), limits = c(0, NA)) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
theme_classic(base_size = 6, base_family = "Arial") +
theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_line(linewidth = 0.2),
legend.position = c(0.91, 0.87),
legend.key = element_blank(),
legend.key.size = unit(1.0, units = 'mm'),
legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
legend.direction = "vertical",
legend.background = element_blank(),
legend.box.background = element_blank(),
legend.box.spacing = unit(-3, 'mm'),
axis.text = element_text(colour = 'black'),
axis.ticks.length = unit(1, units = 'mm'),
axis.ticks = element_line(colour = 'black'),
axis.title = element_text(size = 5),
axis.line = element_line(linewidth = 0.2),
plot.background = element_blank() ) -> p_Dex4
p_Dex4Save to pdf.
ggsave(path = pdf_dir_fig4, filename = "Fig4C_Dex4_vs_WT_Volcano_w_labels.pdf",
plot = p_Dex4, device = cairo_pdf, units = "cm",
height = 5.4, width = 7.0)Prepare data of CSex4 vs WT
res_CSex4 <- subset(res_all, contrast == "CSexon4" & !is.na(p.val))
res_CSex4 <- res_CSex4 |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) Note: nothing is differentially expressed in this contrast
table(res_CSex4$Direction)
None
159
res_KO <- subset(res_all, contrast == "KO" & !is.na(p.val))
res_KO <- res_KO |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) Plot KO vs WT
ggplot(res_KO ) +
aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
geom_point(shape = 21, stroke = 0.2) +
annotate(geom = "label", x = 8.5, y = 0.10, label = "KO",
colour = "white", fill = "gray16", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
annotate(geom = "label", x = -12.5, y = 0.10, label = "WT",
colour = "white", fill = "royalblue3", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
# points on the right side
geom_label_repel(data = subset(res_KO, p.adj <= 0.15 & ratio > 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = -0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
# points on the left side
geom_label_repel(data = subset(res_KO, p.adj <= 0.15 & ratio < 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = 0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
guides(size = guide_legend(override.aes = list(fill = "white"))) +
scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-15, 10), n.breaks = 7) +
scale_y_continuous(expand = expansion(add = c(0, 0.25), mult = 0), limits = c(0, NA)) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
theme_classic(base_size = 6, base_family = "Arial") +
theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_line(linewidth = 0.2),
legend.position = c(0.91, 0.67),
legend.key = element_blank(),
legend.key.size = unit(1.0, units = 'mm'),
legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
legend.direction = "vertical",
legend.background = element_blank(),
legend.box.background = element_blank(),
legend.box.spacing = unit(-3, 'mm'),
axis.text = element_text(colour = 'black'),
axis.ticks.length = unit(1, units = 'mm'),
axis.ticks = element_line(colour = 'black'),
axis.title = element_text(size = 5),
axis.line = element_line(linewidth = 0.2),
plot.background = element_blank() ) -> p_KO
p_KONote: these WT and KO samples were processed at different moment, so there could be some batch effect difference that limits the statistical power.
Save to pdf.
ggsave(path = pdf_dir_fig4, filename = "FigExtra_KO_vs_WT_Volcano_w_labels.pdf",
plot = p_KO, device = cairo_pdf, units = "cm",
height = 12, width = 10)Test against Long Suz12 isoform KO rescue.
data_diff_long <- test_diff(data_norm, type = "control", control = "KO_L")Use same significance thresholds as before.
dep_rescues <- add_rejections(data_diff_long, alpha = 0.05, lfc = 1.5 )Get results to make a volcano plot of the long vs the short
data_results <- get_results(dep_rescues)
data_results |>
as_tibble() |>
select(-significant) |>
pivot_longer(cols = !c("name", "ID"),
names_to = c("contrast", "Feature"),
names_sep = c("_(vs_KO_L_|centered)"),
values_to = "Value"
) |>
mutate(Feature = ifelse(test = Feature == "", yes = "centered", no = Feature)) |>
pivot_wider(names_from = "Feature", values_from = "Value") |>
dplyr::select(-ID) -> res_all_tidy
colnames(res_all_tidy)[colnames(res_all_tidy) == "name"] <- 'First_Col'Add all the info the the tidy results dataframe.
res_all <- left_join(x = res_all_tidy, y = histone_ptms_anno, by = join_by(First_Col))
res_all <- left_join(res_all, WT_area, by = join_by(First_Col) )
res_all |>
mutate(Histone = factor(Histone, levels = sort(unique(res_all$Histone)) ) ) |>
arrange(desc(Histone), Peptide_Start, Peptide_End) |>
ungroup() -> res_allPrepare data for plot.
res_KOrS <- subset(res_all, contrast == "KO_S" & !is.na(p.val))
res_KOrS <- res_KOrS |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) Plot rescues.
ggplot(res_KOrS) +
aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
geom_point(shape = 21, stroke = 0.2) +
annotate(geom = "label", x = 7.5, y = 0.10, label = "KO+S",
colour = "black", fill = "darkorange1", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
annotate(geom = "label", x = -4.85, y = 0.10, label = "KO+L",
colour = "white", fill = "mediumpurple3", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
# points on the right side
geom_label_repel(data = subset(res_KOrS, p.adj < 0.1 & ratio > 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = -0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
# points on the left side
geom_label_repel(data = subset(res_KOrS, p.adj < 0.1 & ratio < 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = 0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
guides(size = guide_legend(override.aes = list(fill = "white"))) +
scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-5.5, 8.25), n.breaks = 7) +
scale_y_continuous(expand = expansion(add = c(0, 0.05), mult = 0), limits = c(0, NA)) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
theme_classic(base_size = 6, base_family = "Arial") +
theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_line(linewidth = 0.2),
legend.position = c(0.10, 0.83),
legend.key = element_blank(),
legend.key.size = unit(1.0, units = 'mm'),
legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
legend.direction = "vertical",
legend.background = element_blank(),
legend.box.background = element_blank(),
legend.box.spacing = unit(-3, 'mm'),
axis.text = element_text(colour = 'black'),
axis.ticks.length = unit(1, units = 'mm'),
axis.ticks = element_line(colour = 'black'),
axis.title = element_text(size = 5),
axis.line = element_line(linewidth = 0.2),
plot.background = element_blank() ) -> p_rescues
p_rescuesSave to pdf.
ggsave(path = pdf_dir_fig4, filename = "FigExtra_KOrS_vs_KOrL_Volcano_w_labels.pdf",
plot = p_rescues, device = cairo_pdf, units = "cm",
height = 6, width = 7.9)I could still add:
rescue efficiency
heatmap with all PTMs.
Now plot to the supplementary figures
There’s no H3.3 signal in the main figure, here I plot only H3.3 signal.
data_K27_36 |>
group_by(Condition, K27, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3) |>
subset(Histone %in% "H3.3") |>
mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) ) |>
ggplot() +
aes(x = Condition, stratum = K27, alluvium = K27,
y = Ratio_Sum, fill = K27) +
geom_alluvium(alpha = 0.5, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, colour = "black") +
geom_fit_text(aes(label = K27 ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = c("K27ac" = "#feac81", "K27me1" ="#ced1af",
"K27me2" = "#748f46", "K27me3" = "#47632a",
"unmod" = "gray" ),
name = 'H3.3') +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm -> p_Alluvial_H33K27
p_Alluvial_H33K27 Save to pdf figure S4D.
ggsave(filename = "FigS4D_H3.3_only_K27_Mrgd_Alluvial.pdf", path = pdf_dir_fig4,
plot = p_Alluvial_H33K27, device = cairo_pdf, units = "cm",
width = 6.5, height = 4.85)Plot the pileup of H3 & H3.3 K36 methylation.
data_K27_36 |>
group_by(Condition, K36, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3) |>
subset(Histone == c("H3", "H3.3") ) |>
mutate(K36 = factor(K36, levels = c("K36me3", "K36me2", "K36me1", "unmod") ) ) |>
ggplot() +
aes(x = Condition, stratum = K36, alluvium = K36,
y = Ratio_Sum, fill = K36) +
facet_wrap(~ Histone, scales = 'free_y') +
geom_alluvium(alpha = 0.5, width = 0.4) +
geom_stratum(alpha = 1, width = 0.6, linewidth = 0.1, colour = "black") +
geom_fit_text(aes(label = K36 ), stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"), size = 7, family = "Arial",
padding.y = grid::unit(0.2, "mm"),
show.legend = F) +
scale_fill_manual(values = c("K36me1" = "#87B7AA", "K36me2" = "#8FA9C2",
"K36me3" = "#8386B7", "unmod" = "gray" ) ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO', 'KO+L', 'KO+S')
) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.title = element_blank()) -> p_Alluvial_K36
p_Alluvial_K36Save to pdf figure S4E
ggsave(path = pdf_dir_fig4, filename = "FigS4E_H3_H3.3_K36_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K36, device = cairo_pdf, units = "cm",
width = 12.0, height = 5)Here I plot extra histone PTMs that were not included in the publication in case someone is interested.
Plot H3K4 modifications.
data_dia_auto |>
group_by(Condition, Modification) |>
subset(Peptide_Start >= 3 & Peptide_End <= 8) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
ggplot() +
aes(x = Condition, stratum = Modification, alluvium = Modification,
y = Ratio_Sum, fill = Modification) +
geom_alluvium(alpha = 0.4, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, size = 0.1) +
geom_fit_text(aes(label = Modification ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = met.brewer(name = "Demuth", direction = 1, n = 5), name = 'H3' ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.text = element_text(margin = margin(l = -9, unit = "mm")),
legend.key.width = unit(6.5, "mm") ) -> p_Alluvial_K4
p_Alluvial_K4Save to pdf.
ggsave(path = pdf_dir_fig4, filename = "FigExtra_H3K4_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K4, device = cairo_pdf, units = "cm",
width = 6.5, height = 4.85)Plot all combinatorial K9 and K14 modifications on histone H3.
data_dia_auto |>
group_by(Condition, Modification) |>
subset(Histone == "H3") |>
subset(Peptide_Start >= 9 & Peptide_End <= 17) |>
mutate(Average_Ratio = mean(Ratio) ) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
ggplot() +
aes(x = Condition, stratum = Modification, alluvium = Modification,
y = Ratio_Sum, fill = Modification) +
geom_alluvium(alpha = 0.4, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, size = 0.1) +
geom_fit_text(aes(label = Modification ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = met.brewer(name = "Hokusai1", direction = -1, n = 15), name = 'H3' ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.text = element_text(margin = margin(l = -12.0, unit = "mm")),
legend.key.width = unit(11.2, "mm"),
legend.key.height = unit(1, 'mm'),
legend.spacing.x = unit(1, 'mm') ) -> p_Alluvial_K9
p_Alluvial_K9ggsave(path = pdf_dir_fig4, filename = "FigExtra_H3K9_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K9, device = cairo_pdf, units = "cm",
width = 6.4, height = 4.85)Plot all combinatorial K27 and K36 modifications on histone H3 and H3.3 variant.
data_dia_auto |>
group_by(Condition, Modification, Histone) |>
subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
ggplot() +
aes(x = Condition, stratum = Modification, alluvium = Modification,
y = Ratio_Sum, fill = Modification) +
facet_wrap(~ Histone) +
geom_alluvium(alpha = 0.4, width = 0.4) +
geom_stratum(alpha = 1, width = 0.6, size = 0.1) +
geom_fit_text(aes(label = Modification ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = met.brewer(name = "Hiroshige", direction = 1, n = 15) ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S'), name = 'H3 & H3.3') +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.text = element_text(margin = margin(l = -14.5, unit = "mm")),
legend.key.width = unit(13.75, "mm"),
legend.key.height = unit(1, 'mm'),
legend.spacing.x = unit(1, 'mm') ) -> p_Alluvial_K27_K36
p_Alluvial_K27_K36This is the best representative picture of K27-K40 modifications relative abundance. This also confirms how the H3.3K36me3 is really the most dramatic change between KO+L and KO+S.
ggsave(path = pdf_dir_fig4, filename = "FigExtra_H3_H3.3_K27_K36_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K27_K36, device = cairo_pdf, units = "cm",
width = 14, height = 5)End analysis.
Protein_PSI_3reps_LS <- read_excel(path = protein_psi_path, sheet = "Sheet1") |>
pivot_longer(cols = starts_with('SUZ12'), names_to = 'Proteoform', values_to = 'PSI') |>
mutate(PSI = PSI * 100)
Protein_PSI_3reps_LS |> group_by(Proteoform) |> summarise(median(PSI))# A tibble: 2 × 2
Proteoform `median(PSI)`
<chr> <dbl>
1 SUZ12-L 90.6
2 SUZ12-S 9.43
Plot
ggplot(Protein_PSI_3reps_LS) +
aes(x = Proteoform, y = PSI, fill = Proteoform) +
geom_boxplot(show.legend = F, outlier.shape = NA, width = 0.5, lwd = 0.1) +
geom_point(shape = 21, size = 0.85, stroke = 0.1, show.legend = F) +
coord_cartesian(ylim = c(0, 100), clip = 'off') +
scale_fill_manual(values = c("SUZ12-L" = "mediumpurple3", "SUZ12-S" = "darkorange1")) +
scale_y_continuous(expand = expansion(mult = 0, add = 0),
n.breaks = 10) +
labs(x = 'Proteoform',
y = 'Rel. abundance in WT ESCs') +
alluvial_thm +
theme(panel.grid.major.y = element_line(colour = 'gray84', linewidth = 0.2),
axis.title.x = element_text(size = 5) ) -> p_prot_PSI
p_prot_PSISave to pdf figure S4C
ggsave(filename = "FigS4C_SUZ12_WT_proteoforms_quantification_barplot.pdf",
plot = p_prot_PSI, device = cairo_pdf, path = pdf_dir_fig4, units = 'cm',
width = 2.65, height = 2.75)sessioninfo::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.1.2 (2021-11-01)
os macOS Catalina 10.15.7
system x86_64, darwin17.0
ui X11
language (EN)
collate C
ctype UTF-8
tz Europe/Madrid
date 2024-02-07
pandoc 2.10.1 @ /usr/local/bin/ (via rmarkdown)
─ Packages ───────────────────────────────────────────────────────────────────
! package * version date (UTC) lib source
ade4 1.7-22 2023-02-06 [1] CRAN (R 4.1.2)
affy 1.72.0 2021-10-26 [1] Bioconductor
affyio 1.64.0 2021-10-26 [1] Bioconductor
annotate 1.72.0 2021-10-26 [1] Bioconductor
AnnotationDbi 1.56.2 2021-11-09 [1] Bioconductor
assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.1.0)
Biobase 2.54.0 2021-10-26 [1] Bioconductor
BiocFileCache 2.2.1 2022-01-23 [1] Bioconductor
BiocGenerics 0.40.0 2021-10-26 [1] Bioconductor
BiocManager 1.30.22 2023-08-08 [1] CRAN (R 4.1.2)
BiocParallel 1.28.3 2021-12-09 [1] Bioconductor
biomaRt 2.50.3 2022-02-06 [1] Bioconductor
Biostrings 2.62.0 2021-10-26 [1] Bioconductor
bit 4.0.5 2022-11-15 [1] CRAN (R 4.1.2)
bit64 4.0.5 2020-08-30 [1] CRAN (R 4.1.0)
bitops 1.0-7 2021-04-24 [1] CRAN (R 4.1.0)
blob 1.2.4 2023-03-17 [1] CRAN (R 4.1.2)
bslib 0.5.1 2023-08-11 [1] CRAN (R 4.1.2)
cachem 1.0.8 2023-05-01 [1] CRAN (R 4.1.2)
Cairo * 1.6-1 2023-08-18 [1] CRAN (R 4.1.2)
callr 3.7.3 2022-11-02 [1] CRAN (R 4.1.2)
cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.1.0)
circlize 0.4.15 2022-05-10 [1] CRAN (R 4.1.2)
cli 3.6.1 2023-03-23 [1] CRAN (R 4.1.2)
clue 0.3-64 2023-01-31 [1] CRAN (R 4.1.2)
cluster 2.1.4 2022-08-22 [1] CRAN (R 4.1.2)
codetools 0.2-19 2023-02-01 [1] CRAN (R 4.1.2)
colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.1.2)
ComplexHeatmap 2.10.0 2021-10-26 [1] Bioconductor
crayon 1.5.2 2022-09-29 [1] CRAN (R 4.1.2)
crosstalk 1.2.0 2021-11-04 [1] CRAN (R 4.1.0)
csaw 1.28.0 2021-10-26 [1] Bioconductor
curl 5.2.0 2023-12-08 [1] CRAN (R 4.1.2)
DBI 1.1.3 2022-06-18 [1] CRAN (R 4.1.2)
dbplyr 2.3.3 2023-07-07 [1] CRAN (R 4.1.2)
DelayedArray 0.20.0 2021-10-26 [1] Bioconductor
DEP * 1.16.0 2021-10-26 [1] Bioconductor
desc 1.4.2 2022-09-08 [1] CRAN (R 4.1.2)
DESeq2 1.34.0 2021-10-26 [1] Bioconductor
devtools 2.4.5 2022-10-11 [1] CRAN (R 4.1.2)
digest 0.6.33 2023-07-07 [1] CRAN (R 4.1.2)
doParallel 1.0.17 2022-02-07 [1] CRAN (R 4.1.2)
dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.1.2)
DT * 0.29 2023-08-29 [1] CRAN (R 4.1.2)
edgeR 3.36.0 2021-10-26 [1] Bioconductor
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.1.0)
evaluate 0.21 2023-05-05 [1] CRAN (R 4.1.2)
fansi 1.0.4 2023-01-22 [1] CRAN (R 4.1.2)
farver 2.1.1 2022-07-06 [1] CRAN (R 4.1.2)
fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.1.2)
fdrtool 1.2.17 2021-11-13 [1] CRAN (R 4.1.0)
filelock 1.0.2 2018-10-05 [1] CRAN (R 4.1.0)
forcats 1.0.0 2023-01-29 [1] CRAN (R 4.1.2)
foreach 1.5.2 2022-02-02 [1] CRAN (R 4.1.2)
foreign 0.8-84 2022-12-06 [1] CRAN (R 4.1.2)
fs 1.6.3 2023-07-20 [1] CRAN (R 4.1.2)
genefilter 1.76.0 2021-10-26 [1] Bioconductor
geneplotter 1.72.0 2021-10-26 [1] Bioconductor
generics 0.1.3 2022-07-05 [1] CRAN (R 4.1.2)
GenomeInfoDb 1.30.1 2022-01-30 [1] Bioconductor
GenomeInfoDbData 1.2.7 2023-08-29 [1] Bioconductor
GenomicRanges 1.46.1 2021-11-18 [1] Bioconductor
GetoptLong 1.0.5 2020-12-15 [1] CRAN (R 4.1.0)
ggalluvial 0.12.5 2023-02-22 [1] CRAN (R 4.1.2)
ggfittext 0.10.0 2023-04-04 [1] CRAN (R 4.1.2)
ggplot2 3.4.3 2023-08-14 [1] CRAN (R 4.1.2)
ggrepel 0.9.3 2023-02-03 [1] CRAN (R 4.1.2)
ggseqlogo 0.1 2017-07-25 [1] CRAN (R 4.1.0)
GlobalOptions 0.1.2 2020-06-10 [1] CRAN (R 4.1.0)
glue 1.6.2 2022-02-24 [1] CRAN (R 4.1.2)
gmm 1.8 2023-06-06 [1] CRAN (R 4.1.2)
gtable 0.3.4 2023-08-21 [1] CRAN (R 4.1.2)
hexbin 1.28.3 2023-03-21 [1] CRAN (R 4.1.2)
hms 1.1.3 2023-03-21 [1] CRAN (R 4.1.2)
htmltools 0.5.6 2023-08-10 [1] CRAN (R 4.1.2)
htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.1.2)
httpuv 1.6.11 2023-05-11 [1] CRAN (R 4.1.2)
httr 1.4.7 2023-08-15 [1] CRAN (R 4.1.2)
impute 1.68.0 2021-10-26 [1] Bioconductor
imputeLCMD 2.1 2022-06-10 [1] CRAN (R 4.1.2)
IRanges 2.28.0 2021-10-26 [1] Bioconductor
iterators 1.0.14 2022-02-05 [1] CRAN (R 4.1.2)
jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.1.0)
jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.1.2)
KEGGREST 1.34.0 2021-10-26 [1] Bioconductor
knitr 1.43 2023-05-25 [1] CRAN (R 4.1.2)
labeling 0.4.2 2020-10-20 [1] CRAN (R 4.1.0)
later 1.3.1 2023-05-02 [1] CRAN (R 4.1.2)
lattice 0.21-8 2023-04-05 [1] CRAN (R 4.1.2)
lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.1.2)
limma 3.50.3 2022-04-07 [1] Bioconductor
locfit 1.5-9.8 2023-06-11 [1] CRAN (R 4.1.2)
magick 2.7.5 2023-08-07 [1] CRAN (R 4.1.2)
magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.1.2)
MALDIquant 1.22.1 2023-03-20 [1] CRAN (R 4.1.2)
MASS 7.3-60 2023-05-04 [1] CRAN (R 4.1.2)
Matrix 1.6-1 2023-08-14 [1] CRAN (R 4.1.2)
MatrixGenerics 1.6.0 2021-10-26 [1] Bioconductor
matrixStats 1.0.0 2023-06-02 [1] CRAN (R 4.1.2)
memoise 2.0.1 2021-11-26 [1] CRAN (R 4.1.0)
metapod 1.2.0 2021-10-26 [1] Bioconductor
MetBrewer 0.2.0 2022-03-21 [1] CRAN (R 4.1.2)
mime 0.12 2021-09-28 [1] CRAN (R 4.1.0)
miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.1.0)
msa 1.26.0 2021-10-26 [1] Bioconductor
MsCoreUtils 1.6.2 2022-02-24 [1] Bioconductor
MSnbase 2.20.4 2022-01-16 [1] Bioconductor
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.1.0)
mvtnorm 1.2-3 2023-08-25 [1] CRAN (R 4.1.2)
mzID 1.32.0 2021-10-26 [1] Bioconductor
mzR 2.28.0 2021-10-27 [1] Bioconductor
ncdf4 1.21 2023-01-07 [1] CRAN (R 4.1.2)
R niar * 0.3.0 <NA> [?] <NA>
norm 1.0-11.1 2023-06-18 [1] CRAN (R 4.1.2)
patchwork 1.1.3 2023-08-14 [1] CRAN (R 4.1.2)
pcaMethods 1.86.0 2021-10-26 [1] Bioconductor
pheatmap * 1.0.12 2019-01-04 [1] CRAN (R 4.1.0)
pillar 1.9.0 2023-03-22 [1] CRAN (R 4.1.2)
pkgbuild 1.4.2 2023-06-26 [1] CRAN (R 4.1.2)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.1.0)
pkgload 1.3.2.1 2023-07-08 [1] CRAN (R 4.1.2)
plyr 1.8.8 2022-11-11 [1] CRAN (R 4.1.2)
png 0.1-8 2022-11-29 [1] CRAN (R 4.1.2)
preprocessCore 1.56.0 2021-10-26 [1] Bioconductor
prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.1.0)
processx 3.8.2 2023-06-30 [1] CRAN (R 4.1.2)
profvis 0.3.8 2023-05-02 [1] CRAN (R 4.1.2)
progress 1.2.2 2019-05-16 [1] CRAN (R 4.1.0)
promises 1.2.1 2023-08-10 [1] CRAN (R 4.1.2)
ProtGenerics 1.26.0 2021-10-26 [1] Bioconductor
ps 1.7.5 2023-04-18 [1] CRAN (R 4.1.2)
psychTools 2.3.6 2023-06-18 [1] CRAN (R 4.1.2)
purrr 1.0.2 2023-08-10 [1] CRAN (R 4.1.2)
R6 2.5.1 2021-08-19 [1] CRAN (R 4.1.0)
rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.1.0)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.1.2)
Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.1.2)
RCurl 1.98-1.12 2023-03-27 [1] CRAN (R 4.1.2)
readr * 2.1.4 2023-02-10 [1] CRAN (R 4.1.2)
readxl * 1.4.3 2023-07-06 [1] CRAN (R 4.1.2)
remotes 2.4.2.1 2023-07-18 [1] CRAN (R 4.1.2)
rjson 0.2.21 2022-01-09 [1] CRAN (R 4.1.2)
rlang 1.1.1 2023-04-28 [1] CRAN (R 4.1.2)
rmarkdown 2.24 2023-08-14 [1] CRAN (R 4.1.2)
rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.1.2)
Rsamtools 2.10.0 2021-10-26 [1] Bioconductor
RSQLite 2.3.1 2023-04-03 [1] CRAN (R 4.1.2)
rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.1.2)
S4Vectors 0.32.4 2022-03-29 [1] Bioconductor
sandwich 3.0-2 2022-06-15 [1] CRAN (R 4.1.2)
sass 0.4.7 2023-07-15 [1] CRAN (R 4.1.2)
scales 1.2.1 2022-08-20 [1] CRAN (R 4.1.2)
seqinr 4.2-30 2023-04-05 [1] CRAN (R 4.1.2)
sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.1.0)
shape 1.4.6 2021-05-19 [1] CRAN (R 4.1.0)
shiny 1.7.5 2023-08-12 [1] CRAN (R 4.1.2)
shinydashboard 0.7.2 2021-09-30 [1] CRAN (R 4.1.0)
stringi 1.7.12 2023-01-11 [1] CRAN (R 4.1.2)
stringr 1.5.0 2022-12-02 [1] CRAN (R 4.1.2)
SummarizedExperiment 1.24.0 2021-10-26 [1] Bioconductor
survival 3.5-7 2023-08-14 [1] CRAN (R 4.1.2)
tibble 3.2.1 2023-03-20 [1] CRAN (R 4.1.2)
tidyr 1.3.0 2023-01-24 [1] CRAN (R 4.1.2)
tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.1.2)
tmvtnorm 1.5 2022-03-22 [1] CRAN (R 4.1.2)
tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.1.2)
urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.1.0)
usethis 2.2.2 2023-07-06 [1] CRAN (R 4.1.2)
utf8 1.2.3 2023-01-31 [1] CRAN (R 4.1.2)
vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.1.2)
vroom 1.6.3 2023-04-28 [1] CRAN (R 4.1.2)
vsn 3.62.0 2021-10-26 [1] Bioconductor
withr 2.5.0 2022-03-03 [1] CRAN (R 4.1.2)
xfun 0.40 2023-08-09 [1] CRAN (R 4.1.2)
XICOR 0.4.1 2023-04-21 [1] CRAN (R 4.1.2)
XML 3.99-0.14 2023-03-19 [1] CRAN (R 4.1.2)
xml2 1.3.5 2023-07-06 [1] CRAN (R 4.1.2)
xtable 1.8-4 2019-04-21 [1] CRAN (R 4.1.0)
XVector 0.34.0 2021-10-26 [1] Bioconductor
yaml 2.3.7 2023-01-23 [1] CRAN (R 4.1.2)
zlibbioc 1.40.0 2021-10-26 [1] Bioconductor
zoo 1.8-12 2023-04-13 [1] CRAN (R 4.1.2)
[1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library
R ── Package was removed from disk.
──────────────────────────────────────────────────────────────────────────────